home *** CD-ROM | disk | FTP | other *** search
/ Power Programmierung / Power-Programmierung (Tewi)(1994).iso / magazine / progjour / 1988 / 01 / whetlin / dlinpakl.for next >
Text File  |  1987-09-20  |  23KB  |  783 lines
  1. *    a TIME function for LAHEY F77L
  2. *
  3. *    Author:    M. Steven Baker
  4. *    Date:    September 20, 1986
  5. *
  6.        real function second()
  7.        integer*4 iticks
  8.        call timer(iticks)
  9.        second = float(iticks)/100
  10.        end
  11.  
  12. *$system
  13.       program linpak
  14.  
  15.       double precision aa(200,200),a(201,200),b(200),x(200)
  16.       double precision time(8,6),cray,ops,total,norma,normx
  17.       double precision resid,residn,eps,epslon
  18.       integer ipvt(200)
  19.       lda = 201
  20.       ldaa = 200
  21. c
  22.       n = 100
  23.       cray = .056
  24.       write(*,1)
  25.     1 format(' Please send the results of this run to:'//
  26.      $       ' Jack J. Dongarra'/
  27.      $       ' Mathematics and Computer Science Division'/
  28.      $       ' Argonne National Laboratory'/
  29.      $       ' Argonne, Illinois 60439'//
  30.      $       ' Telephone: 312-972-7246'//
  31.      $       ' ARPAnet: DONGARRA@ANL-MCS'/)
  32.       ops = (2.0d0*n**3)/3.0d0 + 2.0d0*n**2
  33. c
  34.          call matgen(a,lda,n,b,norma)
  35.          t1 = second()
  36.          call dgefa(a,lda,n,ipvt,info)
  37.          time(1,1) = second() - t1
  38.          t1 = second()
  39.          call dgesl(a,lda,n,ipvt,b,0)
  40.          time(1,2) = second() - t1
  41.          total = time(1,1) + time(1,2)
  42. c
  43. c     compute a residual to verify results.
  44. c
  45.          do 10 i = 1,n
  46.             x(i) = b(i)
  47.    10    continue
  48.          call matgen(a,lda,n,b,norma)
  49.          do 20 i = 1,n
  50.             b(i) = -b(i)
  51.    20    continue
  52.          call dmxpy(n,b,n,lda,x,a)
  53.          resid = 0.0
  54.          normx = 0.0
  55.          do 30 i = 1,n
  56.             resid = dmax1( resid, dabs(b(i)) )
  57.             normx = dmax1( normx, dabs(x(i)) )
  58.    30    continue
  59.          eps = epslon(1.0d0)
  60.          residn = resid/( n*norma*normx*eps )
  61.          write(*,40)
  62.    40    format('     norm. resid      resid           machep',
  63.      $          '         x(1)          x(n)')
  64.          write(*,50) residn,resid,eps,x(1),x(n)
  65.    50    format(1p5e16.8)
  66. c
  67.          write(*,60) n
  68.    60    format(//'    times are reported for matrices of order ',i5)
  69.          write(*,70)
  70.    70    format(6x,'dgefa',6x,'dgesl',6x,'total',5x,'mflops',7x,'unit',
  71.      $         6x,'ratio')
  72. c
  73.          time(1,3) = total
  74.          time(1,4) = ops/(1.0d6*total)
  75.          time(1,5) = 2.0d0/time(1,4)
  76.          time(1,6) = total/cray
  77.          write(*,80) lda
  78.    80    format(' times for array with leading dimension of',i4)
  79.          write(*,110) (time(1,i),i=1,6)
  80. c
  81.          call matgen(a,lda,n,b,norma)
  82.          t1 = second()
  83.          call dgefa(a,lda,n,ipvt,info)
  84.          time(2,1) = second() - t1
  85.          t1 = second()
  86.          call dgesl(a,lda,n,ipvt,b,0)
  87.          time(2,2) = second() - t1
  88.          total = time(2,1) + time(2,2)
  89.          time(2,3) = total
  90.          time(2,4) = ops/(1.0d6*total)
  91.          time(2,5) = 2.0d0/time(2,4)
  92.          time(2,6) = total/cray
  93. c
  94.          call matgen(a,lda,n,b,norma)
  95.          t1 = second()
  96.          call dgefa(a,lda,n,ipvt,info)
  97.          time(3,1) = second() - t1
  98.          t1 = second()
  99.          call dgesl(a,lda,n,ipvt,b,0)
  100.          time(3,2) = second() - t1
  101.          total = time(3,1) + time(3,2)
  102.          time(3,3) = total
  103.          time(3,4) = ops/(1.0d6*total)
  104.          time(3,5) = 2.0d0/time(3,4)
  105.          time(3,6) = total/cray
  106. c
  107.          ntimes = 10
  108.          tm2 = 0
  109.          t1 = second()
  110.          do 90 i = 1,ntimes
  111.             tm = second()
  112.             call matgen(a,lda,n,b,norma)
  113.             tm2 = tm2 + second() - tm
  114.             call dgefa(a,lda,n,ipvt,info)
  115.    90    continue
  116.          time(4,1) = (second() - t1 - tm2)/ntimes
  117.          t1 = second()
  118.          do 100 i = 1,ntimes
  119.             call dgesl(a,lda,n,ipvt,b,0)
  120.   100    continue
  121.          time(4,2) = (second() - t1)/ntimes
  122.          total = time(4,1) + time(4,2)
  123.          time(4,3) = total
  124.          time(4,4) = ops/(1.0d6*total)
  125.          time(4,5) = 2.0d0/time(4,4)
  126.          time(4,6) = total/cray
  127. c
  128.          write(*,110) (time(2,i),i=1,6)
  129.          write(*,110) (time(3,i),i=1,6)
  130.          write(*,110) (time(4,i),i=1,6)
  131.   110    format(6(1pe11.3))
  132. c
  133.          call matgen(aa,ldaa,n,b,norma)
  134.          t1 = second()
  135.          call dgefa(aa,ldaa,n,ipvt,info)
  136.          time(5,1) = second() - t1
  137.          t1 = second()
  138.          call dgesl(aa,ldaa,n,ipvt,b,0)
  139.          time(5,2) = second() - t1
  140.          total = time(5,1) + time(5,2)
  141.          time(5,3) = total
  142.          time(5,4) = ops/(1.0d6*total)
  143.          time(5,5) = 2.0d0/time(5,4)
  144.          time(5,6) = total/cray
  145. c
  146.          call matgen(aa,ldaa,n,b,norma)
  147.          t1 = second()
  148.          call dgefa(aa,ldaa,n,ipvt,info)
  149.          time(6,1) = second() - t1
  150.          t1 = second()
  151.          call dgesl(aa,ldaa,n,ipvt,b,0)
  152.          time(6,2) = second() - t1
  153.          total = time(6,1) + time(6,2)
  154.          time(6,3) = total
  155.          time(6,4) = ops/(1.0d6*total)
  156.          time(6,5) = 2.0d0/time(6,4)
  157.          time(6,6) = total/cray
  158. c
  159.          call matgen(aa,ldaa,n,b,norma)
  160.          t1 = second()
  161.          call dgefa(aa,ldaa,n,ipvt,info)
  162.          time(7,1) = second() - t1
  163.          t1 = second()
  164.          call dgesl(aa,ldaa,n,ipvt,b,0)
  165.          time(7,2) = second() - t1
  166.          total = time(7,1) + time(7,2)
  167.          time(7,3) = total
  168.          time(7,4) = ops/(1.0d6*total)
  169.          time(7,5) = 2.0d0/time(7,4)
  170.          time(7,6) = total/cray
  171. c
  172.          ntimes = 10
  173.          tm2 = 0
  174.          t1 = second()
  175.          do 120 i = 1,ntimes
  176.             tm = second()
  177.             call matgen(aa,ldaa,n,b,norma)
  178.             tm2 = tm2 + second() - tm
  179.             call dgefa(aa,ldaa,n,ipvt,info)
  180.   120    continue
  181.          time(8,1) = (second() - t1 - tm2)/ntimes
  182.          t1 = second()
  183.          do 130 i = 1,ntimes
  184.             call dgesl(aa,ldaa,n,ipvt,b,0)
  185.   130    continue
  186.          time(8,2) = (second() - t1)/ntimes
  187.          total = time(8,1) + time(8,2)
  188.          time(8,3) = total
  189.          time(8,4) = ops/(1.0d6*total)
  190.          time(8,5) = 2.0d0/time(8,4)
  191.          time(8,6) = total/cray
  192. c
  193.          write(*,140) ldaa
  194.   140    format(/' times for array with leading dimension of',i4)
  195.          write(*,110) (time(5,i),i=1,6)
  196.          write(*,110) (time(6,i),i=1,6)
  197.          write(*,110) (time(7,i),i=1,6)
  198.          write(*,110) (time(8,i),i=1,6)
  199.       stop
  200.       end
  201.       subroutine matgen(a,lda,n,b,norma)
  202.       double precision a(lda,1),b(1),norma
  203. c
  204.       init = 1325
  205.       norma = 0.0
  206.       do 30 j = 1,n
  207.          do 20 i = 1,n
  208.             init = mod(3125*init,65536)
  209.             a(i,j) = (init - 32768.0)/16384.0
  210.             norma = dmax1(a(i,j), norma)
  211.    20    continue
  212.    30 continue
  213.       do 35 i = 1,n
  214.           b(i) = 0.0
  215.    35 continue
  216.       do 50 j = 1,n
  217.          do 40 i = 1,n
  218.             b(i) = b(i) + a(i,j)
  219.    40    continue
  220.    50 continue
  221.       return
  222.       end
  223.       subroutine dgefa(a,lda,n,ipvt,info)
  224.       integer lda,n,ipvt(1),info
  225.       double precision a(lda,1)
  226. c
  227. c     dgefa factors a double precision matrix by gaussian elimination.
  228. c
  229. c     dgefa is usually called by dgeco, but it can be called
  230. c     directly with a saving in time if  rcond  is not needed.
  231. c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  232. c
  233. c     on entry
  234. c
  235. c        a       double precision(lda, n)
  236. c                the matrix to be factored.
  237. c
  238. c        lda     integer
  239. c                the leading dimension of the array  a .
  240. c
  241. c        n       integer
  242. c                the order of the matrix  a .
  243. c
  244. c     on return
  245. c
  246. c        a       an upper triangular matrix and the multipliers
  247. c                which were used to obtain it.
  248. c                the factorization can be written  a = l*u  where
  249. c                l  is a product of permutation and unit lower
  250. c                triangular matrices and  u  is upper triangular.
  251. c
  252. c        ipvt    integer(n)
  253. c                an integer vector of pivot indices.
  254. c
  255. c        info    integer
  256. c                = 0  normal value.
  257. c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
  258. c                     condition for this subroutine, but it does
  259. c                     indicate that dgesl or dgedi will divide by zero
  260. c                     if called.  use  rcond  in dgeco for a reliable
  261. c                     indication of singularity.
  262. c
  263. c     linpack. this version dated 08/14/78 .
  264. c     cleve moler, university of new mexico, argonne national lab.
  265. c
  266. c     subroutines and functions
  267. c
  268. c     blas daxpy,dscal,idamax
  269. c
  270. c     internal variables
  271. c
  272.       double precision t
  273.       integer idamax,j,k,kp1,l,nm1
  274. c
  275. c
  276. c     gaussian elimination with partial pivoting
  277. c
  278.       info = 0
  279.       nm1 = n - 1
  280.       if (nm1 .lt. 1) go to 70
  281.       do 60 k = 1, nm1
  282.          kp1 = k + 1
  283. c
  284. c        find l = pivot index
  285. c
  286.          l = idamax(n-k+1,a(k,k),1) + k - 1
  287.          ipvt(k) = l
  288. c
  289. c        zero pivot implies this column already triangularized
  290. c
  291.          if (a(l,k) .eq. 0.0d0) go to 40
  292. c
  293. c           interchange if necessary
  294. c
  295.             if (l .eq. k) go to 10
  296.                t = a(l,k)
  297.                a(l,k) = a(k,k)
  298.                a(k,k) = t
  299.    10       continue
  300. c
  301. c           compute multipliers
  302. c
  303.             t = -1.0d0/a(k,k)
  304.             call dscal(n-k,t,a(k+1,k),1)
  305. c
  306. c           row elimination with column indexing
  307. c
  308.             do 30 j = kp1, n
  309.                t = a(l,j)
  310.                if (l .eq. k) go to 20
  311.                   a(l,j) = a(k,j)
  312.                   a(k,j) = t
  313.    20          continue
  314.                call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
  315.    30       continue
  316.          go to 50
  317.    40    continue
  318.             info = k
  319.    50    continue
  320.    60 continue
  321.    70 continue
  322.       ipvt(n) = n
  323.       if (a(n,n) .eq. 0.0d0) info = n
  324.       return
  325.       end
  326.       subroutine dgesl(a,lda,n,ipvt,b,job)
  327.       integer lda,n,ipvt(1),job
  328.       double precision a(lda,1),b(1)
  329. c
  330. c     dgesl solves the double precision system
  331. c     a * x = b  or  trans(a) * x = b
  332. c     using the factors computed by dgeco or dgefa.
  333. c
  334. c     on entry
  335. c
  336. c        a       double precision(lda, n)
  337. c                the output from dgeco or dgefa.
  338. c
  339. c        lda     integer
  340. c                the leading dimension of the array  a .
  341. c
  342. c        n       integer
  343. c                the order of the matrix  a .
  344. c
  345. c        ipvt    integer(n)
  346. c                the pivot vector from dgeco or dgefa.
  347. c
  348. c        b       double precision(n)
  349. c                the right hand side vector.
  350. c
  351. c        job     integer
  352. c                = 0         to solve  a*x = b ,
  353. c                = nonzero   to solve  trans(a)*x = b  where
  354. c                            trans(a)  is the transpose.
  355. c
  356. c     on return
  357. c
  358. c        b       the solution vector  x .
  359. c
  360. c     error condition
  361. c
  362. c        a division by zero will occur if the input factor contains a
  363. c        zero on the diagonal.  technically this indicates singularity
  364. c        but it is often caused by improper arguments or improper
  365. c        setting of lda .  it will not occur if the subroutines are
  366. c        called correctly and if dgeco has set rcond .gt. 0.0
  367. c        or dgefa has set info .eq. 0 .
  368. c
  369. c     to compute  inverse(a) * c  where  c  is a matrix
  370. c     with  p  columns
  371. c           call dgeco(a,lda,n,ipvt,rcond,z)
  372. c           if (rcond is too small) go to ...
  373. c           do 10 j = 1, p
  374. c              call dgesl(a,lda,n,ipvt,c(1,j),0)
  375. c        10 continue
  376. c
  377. c     linpack. this version dated 08/14/78 .
  378. c     cleve moler, university of new mexico, argonne national lab.
  379. c
  380. c     subroutines and functions
  381. c
  382. c     blas daxpy,ddot
  383. c
  384. c     internal variables
  385. c
  386.       double precision ddot,t
  387.       integer k,kb,l,nm1
  388. c
  389.       nm1 = n - 1
  390.       if (job .ne. 0) go to 50
  391. c
  392. c        job = 0 , solve  a * x = b
  393. c        first solve  l*y = b
  394. c
  395.          if (nm1 .lt. 1) go to 30
  396.          do 20 k = 1, nm1
  397.             l = ipvt(k)
  398.             t = b(l)
  399.             if (l .eq. k) go to 10
  400.                b(l) = b(k)
  401.                b(k) = t
  402.    10       continue
  403.             call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
  404.    20    continue
  405.    30    continue
  406. c
  407. c        now solve  u*x = y
  408. c
  409.          do 40 kb = 1, n
  410.             k = n + 1 - kb
  411.             b(k) = b(k)/a(k,k)
  412.             t = -b(k)
  413.             call daxpy(k-1,t,a(1,k),1,b(1),1)
  414.    40    continue
  415.       go to 100
  416.    50 continue
  417. c
  418. c        job = nonzero, solve  trans(a) * x = b
  419. c        first solve  trans(u)*y = b
  420. c
  421.          do 60 k = 1, n
  422.             t = ddot(k-1,a(1,k),1,b(1),1)
  423.             b(k) = (b(k) - t)/a(k,k)
  424.    60    continue
  425. c
  426. c        now solve trans(l)*x = y
  427. c
  428.          if (nm1 .lt. 1) go to 90
  429.          do 80 kb = 1, nm1
  430.             k = n - kb
  431.             b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
  432.             l = ipvt(k)
  433.             if (l .eq. k) go to 70
  434.                t = b(l)
  435.                b(l) = b(k)
  436.                b(k) = t
  437.    70       continue
  438.    80    continue
  439.    90    continue
  440.   100 continue
  441.       return
  442.       end
  443.       subroutine daxpy(n,da,dx,incx,dy,incy)
  444. c
  445. c     constant times a vector plus a vector.
  446. c     uses unrolled loops for increments equal to one.
  447. c     jack dongarra, linpack, 3/11/78.
  448. c
  449.       double precision dx(1),dy(1),da
  450.       integer i,incx,incy,ix,iy,m,mp1,n
  451. c
  452.       if(n.le.0)return
  453.       if (da .eq. 0.0d0) return
  454.       if(incx.eq.1.and.incy.eq.1)go to 20
  455. c
  456. c        code for unequal increments or equal increments
  457. c          not equal to 1
  458. c
  459.       ix = 1
  460.       iy = 1
  461.       if(incx.lt.0)ix = (-n+1)*incx + 1
  462.       if(incy.lt.0)iy = (-n+1)*incy + 1
  463.       do 10 i = 1,n
  464.         dy(iy) = dy(iy) + da*dx(ix)
  465.         ix = ix + incx
  466.         iy = iy + incy
  467.    10 continue
  468.       return
  469. c
  470. c        code for both increments equal to 1
  471. c
  472. c
  473. c        clean-up loop
  474. c
  475.    20 m = mod(n,4)
  476.       if( m .eq. 0 ) go to 40
  477.       do 30 i = 1,m
  478.         dy(i) = dy(i) + da*dx(i)
  479.    30 continue
  480.       if( n .lt. 4 ) return
  481.    40 mp1 = m + 1
  482.       do 50 i = mp1,n,4
  483.         dy(i) = dy(i) + da*dx(i)
  484.         dy(i + 1) = dy(i + 1) + da*dx(i + 1)
  485.         dy(i + 2) = dy(i + 2) + da*dx(i + 2)
  486.         dy(i + 3) = dy(i + 3) + da*dx(i + 3)
  487.    50 continue
  488.       return
  489.       end
  490.       double precision function ddot(n,dx,incx,dy,incy)
  491. c
  492. c     forms the dot product of two vectors.
  493. c     uses unrolled loops for increments equal to one.
  494. c     jack dongarra, linpack, 3/11/78.
  495. c
  496.       double precision dx(1),dy(1),dtemp
  497.       integer i,incx,incy,ix,iy,m,mp1,n
  498. c
  499.       ddot = 0.0d0
  500.       dtemp = 0.0d0
  501.       if(n.le.0)return
  502.       if(incx.eq.1.and.incy.eq.1)go to 20
  503. c
  504. c        code for unequal increments or equal increments
  505. c          not equal to 1
  506. c
  507.       ix = 1
  508.       iy = 1
  509.       if(incx.lt.0)ix = (-n+1)*incx + 1
  510.       if(incy.lt.0)iy = (-n+1)*incy + 1
  511.       do 10 i = 1,n
  512.         dtemp = dtemp + dx(ix)*dy(iy)
  513.         ix = ix + incx
  514.         iy = iy + incy
  515.    10 continue
  516.       ddot = dtemp
  517.       return
  518. c
  519. c        code for both increments equal to 1
  520. c
  521. c
  522. c        clean-up loop
  523. c
  524.    20 m = mod(n,5)
  525.       if( m .eq. 0 ) go to 40
  526.       do 30 i = 1,m
  527.         dtemp = dtemp + dx(i)*dy(i)
  528.    30 continue
  529.       if( n .lt. 5 ) go to 60
  530.    40 mp1 = m + 1
  531.       do 50 i = mp1,n,5
  532.         dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
  533.      *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
  534.    50 continue
  535.    60 ddot = dtemp
  536.       return
  537.       end
  538.       subroutine  dscal(n,da,dx,incx)
  539. c
  540. c     scales a vector by a constant.
  541. c     uses unrolled loops for increment equal to one.
  542. c     jack dongarra, linpack, 3/11/78.
  543. c
  544.       double precision da,dx(1)
  545.       integer i,incx,m,mp1,n,nincx
  546. c
  547.       if(n.le.0)return
  548.       if(incx.eq.1)go to 20
  549. c
  550. c        code for increment not equal to 1
  551. c
  552.       nincx = n*incx
  553.       do 10 i = 1,nincx,incx
  554.         dx(i) = da*dx(i)
  555.    10 continue
  556.       return
  557. c
  558. c        code for increment equal to 1
  559. c
  560. c
  561. c        clean-up loop
  562. c
  563.    20 m = mod(n,5)
  564.       if( m .eq. 0 ) go to 40
  565.       do 30 i = 1,m
  566.         dx(i) = da*dx(i)
  567.    30 continue
  568.       if( n .lt. 5 ) return
  569.    40 mp1 = m + 1
  570.       do 50 i = mp1,n,5
  571.         dx(i) = da*dx(i)
  572.         dx(i + 1) = da*dx(i + 1)
  573.         dx(i + 2) = da*dx(i + 2)
  574.         dx(i + 3) = da*dx(i + 3)
  575.         dx(i + 4) = da*dx(i + 4)
  576.    50 continue
  577.       return
  578.       end
  579.       integer function idamax(n,dx,incx)
  580. c
  581. c     finds the index of element having max. dabsolute value.
  582. c     jack dongarra, linpack, 3/11/78.
  583. c
  584.       double precision dx(1),dmax
  585.       integer i,incx,ix,n
  586. c
  587.       idamax = 0
  588.       if( n .lt. 1 ) return
  589.       idamax = 1
  590.       if(n.eq.1)return
  591.       if(incx.eq.1)go to 20
  592. c
  593. c        code for increment not equal to 1
  594. c
  595.       ix = 1
  596.       dmax = dabs(dx(1))
  597.       ix = ix + incx
  598.       do 10 i = 2,n
  599.          if(dabs(dx(ix)).le.dmax) go to 5
  600.          idamax = i
  601.          dmax = dabs(dx(ix))
  602.     5    ix = ix + incx
  603.    10 continue
  604.       return
  605. c
  606. c        code for increment equal to 1
  607. c
  608.    20 dmax = dabs(dx(1))
  609.       do 30 i = 2,n
  610.          if(dabs(dx(i)).le.dmax) go to 30
  611.          idamax = i
  612.          dmax = dabs(dx(i))
  613.    30 continue
  614.       return
  615.       end
  616.       double precision function epslon (x)
  617.       double precision x
  618. c
  619. c     estimate unit roundoff in quantities of size x.
  620. c
  621.       double precision a,b,c,eps
  622. c
  623. c     this program should function properly on all systems
  624. c     satisfying the following two assumptions,
  625. c        1.  the base used in representing dfloating point
  626. c            numbers is not a power of three.
  627. c        2.  the quantity  a  in statement 10 is represented to 
  628. c            the accuracy used in dfloating point variables
  629. c            that are stored in memory.
  630. c     the statement number 10 and the go to 10 are intended to
  631. c     force optimizing compilers to generate code satisfying 
  632. c     assumption 2.
  633. c     under these assumptions, it should be true that,
  634. c            a  is not exactly equal to four-thirds,
  635. c            b  has a zero for its last bit or digit,
  636. c            c  is not exactly equal to one,
  637. c            eps  measures the separation of 1.0 from
  638. c                 the next larger dfloating point number.
  639. c     the developers of eispack would appreciate being informed
  640. c     about any systems where these assumptions do not hold.
  641. c
  642. c     *****************************************************************
  643. c     this routine is one of the auxiliary routines used by eispack iii
  644. c     to avoid machine dependencies.
  645. c     *****************************************************************
  646. c
  647. c     this version dated 4/6/83.
  648. c
  649.       a = 4.0d0/3.0d0
  650. C     a = 4.0d0 * 0.33333333d0
  651.    10 b = a - 1.0d0
  652.       c = b + b + b
  653.       eps = dabs(c-1.0d0)
  654.       if (eps .eq. 0.0d0) go to 10
  655.       epslon = eps*dabs(x)
  656.       return
  657.       end
  658.       subroutine mm (a, lda, n1, n3, b, ldb, n2, c, ldc)
  659.       double precision a(lda,*), b(ldb,*), c(ldc,*)
  660. c
  661. c   purpose:
  662. c     multiply matrix b times matrix c and store the result in matrix a.
  663. c
  664. c   parameters:
  665. c
  666. c     a double precision(lda,n3), matrix of n1 rows and n3 columns
  667. c
  668. c     lda integer, leading dimension of array a
  669. c
  670. c     n1 integer, number of rows in matrices a and b
  671. c
  672. c     n3 integer, number of columns in matrices a and c
  673. c
  674. c     b double precision(ldb,n2), matrix of n1 rows and n2 columns
  675. c
  676. c     ldb integer, leading dimension of array b
  677. c
  678. c     n2 integer, number of columns in matrix b, and number of rows in
  679. c         matrix c
  680. c
  681. c     c double precision(ldc,n3), matrix of n2 rows and n3 columns
  682. c
  683. c     ldc integer, leading dimension of array c
  684. c
  685. c ----------------------------------------------------------------------
  686. c
  687.       do 20 j = 1, n3
  688.          do 10 i = 1, n1
  689.             a(i,j) = 0.0
  690.    10    continue
  691.          call dmxpy (n2,a(1,j),n1,ldb,c(1,j),b)
  692.    20 continue
  693. c
  694.       return
  695.       end
  696.       subroutine dmxpy (n1, y, n2, ldm, x, m)
  697.       double precision y(*), x(*), m(ldm,*)
  698. c
  699. c   purpose:
  700. c     multiply matrix m times vector x and add the result to vector y.
  701. c
  702. c   parameters:
  703. c
  704. c     n1 integer, number of elements in vector y, and number of rows in
  705. c         matrix m
  706. c
  707. c     y double precision(n1), vector of length n1 to which is added 
  708. c         the product m*x
  709. c
  710. c     n2 integer, number of elements in vector x, and number of columns
  711. c         in matrix m
  712. c
  713. c     ldm integer, leading dimension of array m
  714. c
  715. c     x double precision(n2), vector of length n2
  716. c
  717. c     m double precision(ldm,n2), matrix of n1 rows and n2 columns
  718. c
  719. c ----------------------------------------------------------------------
  720. c
  721. c   cleanup odd vector
  722. c
  723.       j = mod(n2,2)
  724.       if (j .ge. 1) then
  725.          do 10 i = 1, n1
  726.             y(i) = (y(i)) + x(j)*m(i,j)
  727.    10    continue
  728.       endif
  729. c
  730. c   cleanup odd group of two vectors
  731. c
  732.       j = mod(n2,4)
  733.       if (j .ge. 2) then
  734.          do 20 i = 1, n1
  735.             y(i) = ( (y(i))
  736.      $             + x(j-1)*m(i,j-1)) + x(j)*m(i,j)
  737.    20    continue
  738.       endif
  739. c
  740. c   cleanup odd group of four vectors
  741. c
  742.       j = mod(n2,8)
  743.       if (j .ge. 4) then
  744.          do 30 i = 1, n1
  745.             y(i) = ((( (y(i))
  746.      $             + x(j-3)*m(i,j-3)) + x(j-2)*m(i,j-2))
  747.      $             + x(j-1)*m(i,j-1)) + x(j)  *m(i,j)
  748.    30    continue
  749.       endif
  750. c
  751. c   cleanup odd group of eight vectors
  752. c
  753.       j = mod(n2,16)
  754.       if (j .ge. 8) then
  755.          do 40 i = 1, n1
  756.             y(i) = ((((((( (y(i))
  757.      $             + x(j-7)*m(i,j-7)) + x(j-6)*m(i,j-6))
  758.      $             + x(j-5)*m(i,j-5)) + x(j-4)*m(i,j-4))
  759.      $             + x(j-3)*m(i,j-3)) + x(j-2)*m(i,j-2))
  760.      $             + x(j-1)*m(i,j-1)) + x(j)  *m(i,j)
  761.    40    continue
  762.       endif
  763. c
  764. c   main loop - groups of sixteen vectors
  765. c
  766.       jmin = j+16
  767.       do 60 j = jmin, n2, 16
  768.          do 50 i = 1, n1
  769.             y(i) = ((((((((((((((( (y(i))
  770.      $             + x(j-15)*m(i,j-15)) + x(j-14)*m(i,j-14))
  771.      $             + x(j-13)*m(i,j-13)) + x(j-12)*m(i,j-12))
  772.      $             + x(j-11)*m(i,j-11)) + x(j-10)*m(i,j-10))
  773.      $             + x(j- 9)*m(i,j- 9)) + x(j- 8)*m(i,j- 8))
  774.      $             + x(j- 7)*m(i,j- 7)) + x(j- 6)*m(i,j- 6))
  775.      $             + x(j- 5)*m(i,j- 5)) + x(j- 4)*m(i,j- 4))
  776.      $             + x(j- 3)*m(i,j- 3)) + x(j- 2)*m(i,j- 2))
  777.      $             + x(j- 1)*m(i,j- 1)) + x(j)   *m(i,j)
  778.    50    continue
  779.    60 continue
  780.       return
  781.       end
  782.  
  783.